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ABSTRACT 

A full nova cycle includes mass accretion, thermonuclear runaway resulting in 
outburst and mass loss, and finally, decline. Resumed accretion starts a new cycle, 
leading to another outburst. Multicycle nova evolution models have been calculated 
over the past twenty years, the number being limited by numerical constraints. Here 
we present a long-term evolution code that enables a continuous calculation through 
an unlimited number of nova cycles for an unlimited evolution time, even up to 1.5 x 
10^° yr. Starting with two sets of the three independent nova parameters — the white 
dwarf mass, the temperature of its isothermal core, and the rate of mass transfer on 
to it — we have followed the evolution of two models, with initial masses of 1 Mq and 
0.65 Mq through over 1000 and over 3000 cycles, respectively. The accretion rate was 
assumed constant throughout each calculation: 10~^^ Mq yr~^ for the 1 Mq white 
dwarf, and 10~^ Mq yT~^ for the 0.65 Mq one. The initial temperatures were taken 
to be relatively high: 30 x 10^ K and 50 x 10^ K, respectively, as they are likely to be 
at the onset of the outburst phase. The results show that although on the short-term 
consecutive outbursts are almost identical, on the long-term scale the characteristics 
change. This is mainly due to the changing core temperature, which decreases very 
similarly to that of a cooling white dwarf for a time, but at a slower rate thereafter. 
As the white dwarf's mass continually decreases, since both models lose more mass 
than they accrete, the central pressure decreases accordingly. The outbursts on the 
massive white dwarf change gradually from fast to moderately fast, and the other 
characteristics (velocity, abundance ratios, isotopic ratios) change, too. Very slowly, a 
steady state is reached, where all characteristics, both in quiescence and in outburst, 
remain almost constant. For the less massive white dwarf accreting at a high rate, 
outbursts are similar throughout the evolution. 

Key words: accretion, accretion discs - binaries: close - novae, cataclysmic variables 
- white dwarfs 



1 INTRODUCTION 

It is now twenty years since the publication of the first 
numerical simulation of a full nova cycle (Prialnik, 1986). 
During this time modelling of nova outbursts has advanced 
significantly and calculations have extended to continuous 
evolution thorough a series of several cycles (e.g., Shara et 
al., 1993). Extensive parameter studies of the classical nova 
phenomenon have been carried out by Prialnik and Kovetz 
(1995) and Kovetz and Prialnik (1997), and more recently 
by Yaron et al. (2005). In parallel, computation capability 
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has increased tremendously and thus the computing time re- 
quired to simulate a full nova cycle has dropped from months 
to minutes. 

Models have shown — and observations have confirmed 
— that the characteristics of a nova outburst are deter- 
mined by three independent parameters (Schwartzman et 
al., 1994; Prialnik and Kovetz, 1994; Prialnik, 1995). On 
the theoretical side, these are the mass and core temper- 
ature (or luminosity) of the white dwarf (WD), which ac- 
cretes mass from a binary companion and eventually ignites 
it explosively and ejects it back into space, and the rate of 
accretion on to the WD. Changing these parameters inde- 
pendently over the entire ranges allowed by stellar structure 
theory reproduces the wide range of observed characteristics 
of novae (Yaron et al., loc. cit.). So far, therefore, the theory 
of nova outbursts has proved extremely successful. However, 
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the three free parameters are not truly independent. They 
are hnked by the long-term evolution of the binary system, 
which changes all of them simultaneously. Thus the question 
remains whether the entire 3-dimensional parameter space 
is accessible to nova binary systems. Is it possible for the 
WD to cool to low temperatures despite repeated outbursts 
and loss or gain of mass? Does the WD mass remain the ex- 
clusive structural parameter, or is the structure affected on 
the long-term scale by the events occurring at the surface? 
Such questions can only be answered by following the evo- 
lution of an accreting WD continuously, not through one or 
a few nova cycles, but through a thousand such cycles and 
more. This is the purpose of the present paper. Long-term 
unabridged evolutionary calculations through repeated nova 
outburst cycles require some changes in the evolution code 
used in previous calculations, which will be described in Sec- 
tion (2] The results of two computation runs are described 
and analysed in Section |31 both from the point of view of 
the outbursts and of the WD on top of which they occur. 
A brief discussion, summary and conclusions are given in 
Section g] 



2 LONG-TERM EVOLUTIONARY 
CALCULATIONS 

The hydrodynamic Lagrangian stellar evolution code used 
for the calculations presented here was described in some 
detail by Prialnik and Kovetz (1995). It includes OPAL 
opacities, a nuclear reactions network among 40 heavy ele- 
ment isotopes up to ^'^P, convection according to the mixing- 
length theory, diffusion for all elements, accretional heating, 
and a mass-loss algorithm that applies a steady, optically 
thick supersonic wind solution. We note, in particular, that 
the dynamical phases are calculated by solving the equation 
of motion along with the energy balance equation (rather 
than imposing hydrostatic equilibrium). Mass loss is calcu- 
lated continuously, according to the mass loss rate Mm-i 
derived from the optically thick wind solution. At each time 
step 5t, an amount Mm~i5t is subtracted from the outer- 
most mass shell Ams. Time steps are constrained during 
this phase by the requirement Mm-iSt < Ams- Whenever 
Ams becomes very small, it is merged with the underlying 
mass shell. In a similar manner, during the accretion phase, 
an amount MaccSt is added to the outermost mass shell, 
and whenever Ams exceeds a prescribed value, the shell is 
divided into two mass shells. 

Although thousands of cycles are computed in the 
present study, each nova cycle, including all evolutionary 
phases and involving a large network of nuclear reactions, is 
calculated accurately, with no short-cuts or simplifications. 
A detailed description of the changing physical parameters 
during the evolution of a full nova outburst was given by Pri- 
alnik (1986). Here we refrain from dealing with any single 
outburst cycle in detail, but focus on the trends of change of 
typical characteristics and on long-term evolutionary effects. 

The spherical grid adopted for the numerical solution, 
which includes the entire WD, is fixed, except for the outer- 
most grid shells, which may grow or shrink, and be divided 
or merged, during the accretion and mass-loss episodes of 
the nova cycle. The fixed grid is, however, by no means 
an equally-spaced one. The masses of shells vary by many 



orders of magnitude, according to the gradients of various 
physical properties (such as temperature, density, element 
abundances) that are expected to develop during evolution. 
Thus, the spatial grid must be prepared in advance in an- 
ticipation of the evolutionary course of the stellar model. In 
particular, the outer layers of the WD that are expected to 
take part in the various outburst processes (diffusion, con- 
vection, nuclear burning, acceleration, expansion, ejection, 
contraction, and so forth) must be finely zoned. Since, as 
a rule, each outburst erodes the WD of some mass (in ad- 
dition to the accreted mass), the grid is set in advance for 
a prescribed number of nova cycles. The larger the number 
of cycles, the larger the grid, and consequently, the larger 
the computation time. This was the procedure adopted in 
former multicycle calculations. Thus the number of cycles 
was restricted in the original form of the code by the fixed 
spacial grid of the numerical scheme. 

The long-term evolution code should allow an unlim- 
ited number of cycles. A new algorithm is therefore devised 
whose purpose is to rezone the grid at the beginning of each 
nova cycle, thereby preparing it for the upcoming evolution- 
ary stages of that cycle. This numerical process, which in- 
volves interpolation on a previous grid, is designed in such a 
way as to conserve energy and mass of each nuclear species. 
Hydrostatic equilibrium (which prevails between the dynam- 
ical phases of the outburst itself) is preserved as well. 

One of the problems that have to be overcome is the de- 
termination of the optimal mass shell size for the outer part 
of the grid. Since consecutive nova outbursts are very similar 
(almost identical), this is done at each cycle based on the 
characteristics of the previous one, namely, the amount of ac- 
creted mass and ejected mass, and the depth of the burning 
shell and convective zone. Another problem is to determine 
the depth of the finely zoned region and the rate at which to 
gradually increase shell masses until they eventually match 
the coarsely zoned bulk of the star. This is achieved by a 
geometrically increasing series, with a typical factor of 1.25. 
Since matrix inversion is required at each iteration of every 
time-step of the numerical calculation, a compromise must 
be reached between computing time and resolution (or ac- 
curacy). Typically, the grid included between 200 and 400 
mass shells for the entire configuration and each continuous 
— hands-off — run took several weeks of CPU time on a 
fast PC. 

Finally, the code is capable of correcting itself in case 
of failure, by going back to an earlier stage and continuing 
from that stage with adjusted numerical parameters, such as 
length of time-step, maximal number of iterations allowed 
for achieving convergence, and so forth. Occasionally, a cou- 
ple of times per thousand cycles, it does happen that the 
code goes astray, but it recovers without intervention. Given 
the complexity of the input physics, in particular the need 
to interpolate among a large number of opacity tables, this 
should not be surprising. Such stray points have been re- 
moved from the results. 

A run spanning thousands of nova outburst cycles pro- 
duces a prohibitive amount of data, which must be carefully 
selected, divided among a large number of files, stored, and 
analysed, mainly by graphical means. A great deal of in- 
formation is necessarily lost in the process and can only be 
retrieved by repeating the calculation. 
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Figure 1. Evolution of Model A (MwD = ^Mq, M = 
1O-"M0 yr-i and Two = 30 x 10"^ K): top left: ac- 
creted/ejected/envelope mass; top right: maximal temperature 
attained at outburst; bottom left: ejected/envelope metallicity; 
bottom right: ejected/envelope helium mass fraction. 



3 EVOLUTION THROUGH REPEATED NOVA 
CYCLES 

Two evolution runs were carried out, starting from two dif- 
ferent initial parameter combinations taken from the ex- 
tended grid of nova models studied by Yaron et al. (2005). 
We chose a relatively massive WD, of relatively high tem- 
perature, accreting at a low rate (hereafter. Model A), and a 
hot, low-mass WD, accreting at a relatively high rate (here- 
after. Model B). According to Prialnik and Kovetz (1995) 
and Yaron et al. (2005), the outburst characteristics of these 
models differ considerably, although in both cases they were 
found to be well within the observed range of nova proper- 
ties. 
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Figure 2. Evolution of Model A [MwD = 1 Mq ,M = 
10"" Mq yr-i and Two = 30 x 10'^ K): top: average veloc- 
ity of the ejecta Vavg] bottom: duration of the mass loss phase 



3.1 Nova outbursts — Model A 



The initial parameters of Model A are: Mwd ~ 1 Mq, M = 
10"" Mq yr"^ and Two = 30 x 10® K. Evolution was 
followed for 1.5 x 10^" yr (roughly a Hubble time), during 
which time the WD underwent 1001 nova outbursts. 

The long-term evolution of the main physical charac- 
teristics related to the outbursts that occur on top this ac- 
creting WD is shown in Figures [T]|31 The top left panel of 
Fig. [T]shows the accreted mass mace, ejected mass m^j and 
the mass of the convective envelope that develops after hy- 
drogen ignition m^nv . All three increase with time for the 
first ~240 cycles, or ~ 10^ yr, but stabilize thereafter to a 
very slow rate of growth. At first, the ejected mass is al- 
most identical to the envelope mass, which means that the 
remnant, hydrogen-rich layer at the end of the outburst is 
very small, and thus the time of decline of the bolometric 
magnitude is short as well. We recall that the decline starts 
when the remnant hydrogen is exhausted. Gradually, how- 
ever, the trend shifts towards the envelope mass being al- 
most equal to the accreted mass, and the ejected mass being 
larger than both. This means that WD material is ejected 
directly during the late stages of mass loss, without being 



previously mixed with the accreted material. Observation- 
ally, this should be reflected by a change in composition of 
the ejecta, which may be detectable if/when the material 
ejected earlier becomes transparent. The reason for this ef- 
fect is a change within the evolution of the nova outburst. 
When most of the mass contained in the convective enve- 
lope has been ejected, the star contracts back to almost 
WD size. The high (bolometric) luminosity persists for a 
while longer, until the hydrogen in the envelope remnant 
is burnt out. When the remnant is massive, however, con- 
traction followed by heating of the hydrogen-rich remnant 
layer may ignite it explosively again, which will lead to a 
second mass loss episode, so that the entire original enve- 
lope and even some additional WD material will be ejected. 
Such behaviour has already been encountered and discussed 
in previous nova outburst studies (e.g., Prialnik and Livio 
1995). 

For each outburst, we calculate the average abundance 
of any element n in the ejecta Xn^ej by 



Xn,ej = J Mm-lXn,s{t)dt/mej, (1) 

where X„,g(t) is the respective abundance in the outermost 
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mass shell from which mass is removed (see Section |2} at 
a given time; it is not necessarily constant (see also Kovetz 
and Prialnik 1997), as it depends upon the mixing history 
of the envelope following the thermonuclear runaway. The 
integration is carried over the entire mass loss phase of the 
cycle. The ejecta average metallicity Z^j decreases with re- 
peated cycles, as the accreted mass increases. At the same 
time the helium mass fraction Yej increases, but both reach 
a plateau after 240 cycles, when the accreted and ejected 
masses stabilize as well. At this stage, the breakdown of Zej 
into the most abundant elements is: N - 0.09, O - 0.08, and 
C - 0.04. The difference between Z^j and Z^nv is a direct 
consequence of the relationship between niej and nienv, ex- 
plained above: Zej > Zenv, when during the secondary mass 
loss episode, WD material {Z ~ 1) that was not previously 
mixed with the accreted material, is ejected as well. 

The maximal temperature attained at outburst reaches 
a maximum of about 2.3 x 10* K around the 400'th cycle 
and declines very slowly thereafter. 

Perhaps the most interesting result of this long-term 
calculation is the relatively sharp transition from fast nova 
outbursts — for roughly the first 200 cycles — to moder- 
ately fast nova outbursts from the 250 'th outburst onwards, 
as shown by the duration of the mass loss episode plotted 
in the lower panel of Fig. (2] Even more intriguing is the be- 
haviour of the average expansion velocity, shown in the up- 
per panel: as the fast nova outbursts progress, the velocity 
goes through a maximum and settles to an almost constant 
value of about 500 km s~^ for the moderately fast outbursts. 
It thus appears that fast novae may exhibit a range of ex- 
pansion velocities ranging between 500 and 2500 km s~^, 
while slower ones should have typically lower velocities. We 
should bear in mind, however, that a low average velocity 
may still allow much higher peak velocities for short peri- 
ods of time. Of course, these conclusions relate to outbursts 
obtained for a 1 Mq WD and may change with the WD 
mass. 

The isotopic ratios shown in Fig. |3] undergo a similar 
evolution: a changing, non-monotonic pattern along the first 
250 fast nova cycles, settling down to constant ratios for 
most of the evolution. We note significant overabundances 
of ^^C and ^^O, compared to solar abundances (taken from 
Lodders, 2003), a somewhat lower overabundance of ^^N, 
while is underabundant. We also note that for the first 
few tens of outbursts the trends of change are different and 
variable. Thus, in this respect, results based on one or even a 
few nova cycles may be misleading. In other respects, how- 
ever, a general trend is already apparent and maintained 
from the beginning. 

3.2 Nova outbursts - Model B 

The initial parameters of Model B are: Mwd ~ 0.65 Mq, 
M = 10"® Mq yr"^ and Two = 50 x lO'' K. Evolution was 
followed for a much shorter time, ~ 5 x 10* yr, but dur- 
ing this time the WD underwent over 3000 nova outbursts. 
The total mass transferred to it by accretion amounted to 
~ 0.5 Mq , which is of the order of a typical red dwarf binary 
companion mass. Thus, in spite of the relatively short evo- 
lution time, the companion mass may have been exhausted 
by its end. 

The long-term evolution of the main physical character- 
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Figure 3. Evolution of isotopic ratios normalized by the cor- 
responding solar ratios for Model A: [i^C/i^C], [ISn/i^N], 
[i^o/iSO], [i^O/i^o]. 
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Figure 4. Evolution of Model B (MwD = 0.65 Mq, M = 
IO^^Mq yr-i and Two = 50 x lO'^ K): top left: ac- 
creted/ejected/envelope mass; top right: maximal temperature 
attained at outburst; bottom left: ejected/envelope metallicity; 
bottom right: ejected/envelope helium mass fraction. 

istics related to nova outbursts that occur on top this accret- 
ing WD is shown in Figures [3] and [S] Here, due to the high 
accretion rate, the evolutionary timescale is much shorter; 
thus, for example, 1000 cycles which spanned ~ 1.5 x 10^" yr 
for Model A, take only ~ 10* yr for model B. The times 
are not precisely inversely proportional to the correspond- 
ing accretion rates because the accreted masses required for 
outbursts to be triggered depend on the WD mass. 

The striking difference between Models A and B, which 
is immediately obvious from the comparison of Figs. [1] and 
|3]as well as Figs. [3] and [S] is the constancy (or monotony) in 
the evolution of the outburst characteristics for the low-mass 
WD, in contrast to the sharp changes that take place in the 
course of evolution of the more massive one. But, in fact, 
indication for this trend of behaviour was already provided 
by the grid of nova models (Yaron et al. 2005). We note. 
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Figure 5. Evolution of isotopic ratios normalized by the cor- 
responding solar ratios for Model B: [l^C/i^c], [i^N/i^N], 
pO/i«0], [iSO/i^O]. 

in particular, that for Model B the envelope mass is always 
larger than both the accreted and the ejected masses. The 
mass fractions of the most abundant elements of the ejecta 
are 0.065 for N, 0.055 for O, and 0.004 for C. Regarding 
isotopic ratios, we note that ^^N is underabundant in the 
present case, whereas for Model A it was overabundant at 
various levels throughout most of the evolution. The aver- 
age expansion velocity oscillates between 150 km s~^ and 
220 km s~^, with no particular trend. The outburst du- 
ration, as measured by the extent of the mass loss phase, 
increases gradually from ~80 to ^^160 days along the first 
1000 cycles and oscillates between these limits thereafter. 

3.3 Evolution of the WD - Models A and B 

One of the most puzzling questions related to nova outbursts 
is their effect, if any, on the characteristics of the underly- 
ing WD, keeping in mind that the outbursts are confined to 
the outermost layers of the WD, less than a thousandth of 
its mass. Previous studies, even if they considered several 
cycles, did not span a sufficient fraction of the WD evolu- 
tionary timescale for answering this question. An indication 
that the WD core may continue to cool despite the heat gen- 
erated by nova outbursts was given by Prialnik (1987). In 
the present calculations the evolution of the WD was carried 
out for 1.5 X 10^" yr (about a Hubble time); we are thus in 
a position to answer this question. 

Using the same code, we followed the evolution of the 
same initial WD configuration, omitting accretion, to ob- 
tain the cooling curve of the WD. The central temperatures 
of the undisturbed WD and of the accreting one are com- 
pared in Fig. m For the 1 M© WD (Model A), the accreting 
WD cools at the same rate as the single one for a period of 
^ 10^ yr (equivalent to about 250 cycles), but then the cool- 
ing trend slows down and the WD core tends to a constant 
temperature of ~ 5 x 10® K. The evolutionary timescale for 
Model B is short compared to the cooling timescale of a 
WD, hence the cooling curves are still close even after sev- 
eral thousand cycles, which amount to only a few 10* yr. 
The luminosities of the cooling WD and of the accreting 
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Figure 6. Evolution of the central temperature (top) and of the 
luminosity [bottom) of the accreting WD and a cooling WD of 
the same mass and initial core temperature: left — Model A; right 
— Model B. The accretion luminosity, eq. l[2]l, is also shown in 
the lower panels. An additional curve for Model B shows the 
long-term evolution through almost 3000 nova cycles, adopting 
the same mass and accretion rate, but a 10 times lower initial 
temperature, Two = 5 X 10® K (see discussion in Section l4.1ll . 



WD in quiescence (during the accretion episodes separat- 
ing outbursts) are compared in Fig. [S] The luminosity of 
the accreting WD has two sources: the heat emanated by 
the coofing WD core and the energy imparted by accretion, 
known as accretion luminosity, given by 

_ aaccGMwoM 

l^acc — 5 (Zj 

HWD 

where Uacc is taken to be 0.15, following Regev and Shara 
(1989). The outer layers of the WD adjust so as to re-emit 
this influx of energy, in a similar manner to an irradiated 
star (Kovetz et al. 1988). We note that only at the begin- 
ning, while the WD core is still relatively hot, does the core 
luminosity exceed the accretion luminosity. 

In order to understand the thermal evolution of the 
WD, we plot in Fig. [7] temperature profiles at various times 
during evolution (marked by the cycle number). All pro- 
files correspond to the end of a nova cycle, just before the 
onset of a new accretion episode. Thus the configuration 
includes the WD and the eventual remnant of the enve- 
lope at the end of the mass loss phase. At first, for several 
hundred cycles, the temperature gradient is negative every- 
where. Gradually, however, the outer layers of the WD, well 
beneath the outermost one that is directly affected by the 
outburst (mainly through diffusion, which changes its com- 
position), are heated and a temperature inversion develops. 
Nevertheless, for many additional cycles, this temperature 
"wall" does not prevent the core from cooling, since the tem- 
perature still decreases between the centre and the inversion 
zone. But the slope of the core temperature profile slowly de- 
creases and finally flattens. At this stage the WD core ceases 
to cool. It may subsequently start to heat up, but for this 
to happen, more than a Hubble time seems to be required 
for the massive WD, and several thousand additional nova 
cycles (that is, a massive companion) for the less massive 
one. In conclusion, the shape and evolution of the tempera- 
ture profiles explains the difference (or similarity) between 



6 N. Epelstain 0. Yaron A. Kovetz and D. Prialnik 



Time(IO^yr) Time(1o''yr) 




log(Depth/R^^^) time[yr] ^ time[yr] 




Figure 7. Temperature profiles from the surface inwards on a 
logarithmic scale at the end of a cycle (onset of accretion) for 
several cycles — as marked — along the evolutionary course of 
Model A {top) and of Model B (bottom), luminosity, eq. {J}, is 
also shown: Ze/t - Model A; right - Model B. 



the cooling curves of an accreting and an undisturbed WD, 
shown in Fig. |6] 

The evolution of the central density of the accreting WD 
is shown in the lower panel of Fig. [S] together with that of 
the single WD. While the latter increases approaching an 
asymptotic value, the former — for Model A — increases 
at first but decreases steadily thereafter. The reason for this 
behaviour lies in the competition between the already men- 
tioned cooling of the core, which drives the density up, and 
decrease of the WD mass (shown in the upper panels of 
Fig. (Sjl, which drives it down. Applying to the WD centre 
the equation of state for degenerate matter to first order in 
temperature (e.g.. Landau and Lifshitz, 1969), we have 



Pc = Pp. 



5/3 



4/3 



where, in standard notation. 



/3 = 



h'^ (Stt 
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Figure 8. Decrease of the WD mass due to mass loss during nova 
outburst cycles (top) and evolution of the central density of the 
accreting WD and a cooling WD of the same mass and initial 
core [bottom) temperature: left - Model A; right — Model B. 



and fie ~ 2. 

On the other hand, hydrostatic equilibrium for a poly- 
tropic equation of state requires 



(5) 



where B„ is a constant determined by the polytropic index n 
(e.g., Chandrasekhar, 1967), thus B1.5 = 0.206 for P oc p^/^. 
Combining eqs.|3]and[S] we obtain a relation between pc, Tc 
and MwD, of the form 



/3pJ/3 ^ ^ ^^^^ 

pc 

where 

7 = (47r)^''^0.206G 

Taking the time derivative of eq. (O, we obtain 

aT^\dpc ^ndT.-y i/3dMwD 



f2/3 



(6) 



(7) 



2pc 



1/3 

Pc 



Tc dTc -y ,_j 1^ 



Since a is a very small number, deriving from the small 
correction to the degenerate equation of state due to temper- 
ature, the coefficient of the density derivative on the LHS of 
eq. ([H is positive, regardless of the changes with time in the 
temperature and density values. The terms on the RHS have 
opposite signs, since both dT^/dt and dMwo/dt are nega- 
tive. Since Tc changes rapidly with time at the beginning 
and very slowly thereafter, while the WD mass decreases 
at an almost constant rate, the LHS will change sign from 
positive to negative at some point and thus the density will 
go through a maximum. The competition between the two 
terms is shown in Fig. [5] for Model A. A correction factor 
was introduced in the second term, to take account of the 
fact that for a 1 Mq WD relativistic effects on the equa- 
tion of state ((31) are non- negligible. The density maximum 
corresponds to the intersection. 

However, since the rate of change of the WD mass is 
roughly proportional to the rate of accretion, 
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Figure 9. Left: Evolution of tiie LHS terms of eq. JSj for Model A 
(see text for details): Fl represents the first term, which includes 
the derivative of temperature, and F2 - the second term, which 
includes the derivative of the WD mass. Right: Zoom in into the 
evolution of the central density (shown in Fig.[8ll, to show the rise 
and decline of pc during early evolution. 



dMwD 
dt 



-M 



(9) 



for a very high accretion rate the second term on the RHS 
of eq. (|8]l will always dominate and the central density will 
decrease monotonically. This is the situation for Model B, 
where the accretion rate is two orders of magnitude higher 
than for Model A. Equally, if dTc/dt > 0, the central density 
will decrease monotonically. 



3.4 Energy budget 

It is instructive to examine the energy budget of nova out- 
bursts and its eventual change with time. The nuclear energy 
generated during the thermonuclear runaway and during 
the equilibrium burning of the remnant hydrogen (if any), 
Enuc, is spent in radiation at close to (or exceeding) Edding- 
ton luminosity during the outburst, Erad, kinetic energy of 
the ejecta, Ekin, and gravitational energy required to lift 
the ejected material out of the potential well of the WD, 
Egrav These are plotted in Fig. [TD] for Model A. The last 
term is not computed during calculations, but estimated by 
Egrav ~ MwDTn^j / R{Mwd), which is an upper limit to 
the actual value, since R{Mwd) is the radius at the lower 
boundary of the ejected mass. We note that the bulk of nu- 
clear energy is spent in work against the gravitational field 
of the WD. We also note that the radiated energy exceeds 
the kinetic energy and, for most of the evolution, quite sig- 
nificantly. Finally, the sum of sinks is lower than the source, 
that is, 



Enuc ^ Egrav ~\~ Erad ~t~ Ej^ 



(10) 



which means that a small fraction of the heat generated 
at outburst is conducted into the WD. We may estimate 
the total thermal energy that the WD gains for the entire 
evolution time, or more precisely, a lower limit for it, by 

J2aycles ^"^'^ " i^Srav + Erad + Ekir.) > 2 X 10*^ Crg. This 

energy represents about 1% of the total energy of the WD 
and explains the very slow change in slope of the core tem- 
perature profile shown in Fig. [7] During evolution, the WD 
expands due to mass loss (see Fig. IS}, and therefore its grav- 
itational potential (binding) energy — as computed by the 
evolution code — increases (that is, becomes less negative), 
while the internal energy (mostly the energy of the degen- 
erate electron gas) decreases by roughly the same amount 
(~ 15%), as expected (cf. Mestel and Ruderman 1967). 
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Figure 10. Energy sources and sinks as a function of time: total 
nuclear energy released during a nova cycle, gravitational binding 
energy of the ejecta, total radiated energy during outburst and 
kinetic energy of the ejecta, for Model A. 



For the 0.65 M© WD, the budget is slightly different. 
All energy values are almost constant in time, as are most 
of the other nova characteristics. Typical values are: E„uc ~ 
7.6 X Iff^s erg, Egrav ~ 2.7 x lO'*^ erg, Erad « 4.2 x 10^'' erg, 
and finally, Ekin ~ 8 x 10*^ erg. We note that in this case the 
radiated energy is the major energy sink — due to the long 
duration of the outburst — while the total kinetic energy is 
only a very small fraction of the nuclear energy generated, 
considerably lower than the radiated energy. 



4 DISCUSSION AND CONCLUSIONS 

4.1 Comparison with earlier studies 

Model A starts at a core temperature of 30 x 10® K, which 
decreases during the WD's evolution. After 210 cycles (7.8 x 
10** yr), the temperature reaches 10 x 10® K. This result 
enables a comparison with the grid of nova models calculated 
by Prialnik and Kovetz (1995) and Yaron et al. (2005), where 
the different values of WD core temperature were set as 
initial independent parameters. Similarly, the 0.65 Mq WD, 
which started at a core temperature of 50 x 10® K, reaches 
30 X 10® K after 970 cycles (9.7 x lO'' yr). The resuhs are 
compared in Table [T] We note a striking similarity between 
the results, allowing for the fact that the WD mass changes 
as well along with the core temperature. In the case of the 
0.65 Mq WD, for example, it dropped to 0.637 Mq by the 
time the core temperature dropped to 30 x 10® K. 

Townsley and Bildsten (2004) studied the problem of 
accreting WDs analytically, assuming that the WD core 
reaches an equilibrium temperature and this is achieved by 
heating of the core due to release of gravitational energy 
during compression of the accreted material as well as quiet 
nuclear burning during the nova cycle. Their study, by its 
nature, does not take into consideration evolutionary effects, 
such as the decrease of the WD mass, and dynamical phases 
of the outburst, which determine to a large extent the en- 
ergy budget. Our full-scale calculations show that, indeed. 
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Table 1. Compajrison between present calculations and the grid of nova models (Yaron et al., 2005) for M^yd = 1 and M = 
10~^^Mq yr~i and for MwD = 0.65 Mq and M = 10~^Mq yr~^ - Characteristics of the nova outbursts. 
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a steady state (or equilibrium) is eventually attained or ap- 
proached asymptotically, but on a very long timescale. 

The equilibrium core temperatures that we find for the 
two cases considered here are considerably higher than those 
derived by Townsley and Bildsten (loc. cit.); the accreted 
masses, equivalent to their ignition masses, for the same 
cases, are similar. In order to test whether the steady state 
temperature is, indeed, an equilibrium temperature, we have 
run a calculation similar to Model B, but for a very low 
initial WD temperature, lower than the asymptotic value. 
The evolution was again followed for close to 3000 full nova 
cycles and we found the WD temperature to rise slowly to- 
wards the same steady state value of about 1.5 x 10^ K, as 
the cooling, initially hot, WD. We wish to stress, however, 
that in both cases the asymptotic value is approached very 
slowly compared with the evolutionary timescale of the nova 
binary system, on which the mass ratio of the components 
changes very significantly, and the mass transfer rate is ex- 
pected to change as well. Consequently, the equilibrium core 
temperature may not be assumed to occur a priori. 



4.2 Summary 

We have carried out long-term continuous evolutionary cal- 
culations of accreting WDs through thousands of nova cy- 
cles, where each cycle was followed in detail through its dif- 
ferent phases: accretion, thermonuclear runaway, expansion, 
mass loss and contraction. 

Ths study shows, for the first time, the long-term effect 
nova outbursts have on the structure and evolution of the 
WD at the surface of which they occur. At the same time, it 
shows the eflFect of the changing WD characteristics, in turn, 
on the nova outbursts. We have chosen two very different 
combinations of WD mass and accretion rate in order to test 
these two types of effect. The main results and conclusions 
reached may be summarized as follows: 

(i) The WD physical parameters change considerably 
with time in both cases, in particular: 

(a) The WD cools as a single, undisturbed WD for a 
while, but after a time, cooling slows down and will even- 
tually cease. The maximum temperature attained in the 
burning shell Tmax increases at first, but settles at an al- 
most constant value as cooling of the WD slows down. 

(b) The WD luminosity decreases with time, but in 
contrast to a single WD, only down to the accretion lumi- 
nosity value, where it remains constant. This, of course, 
applies only to the quiescent phases of accretion between 
outbursts. 



(c) The WD mass decreases steadily for the cases con- 
sidered, but not by much. The mass of the donor star, by 
contrast, decreases considerably: by 0.16 Mq for Model 
A, and by 0.51 Mq for Model B. 

(d) The WD becomes less dense as its mass decreases, 
as expected from the relation pc{M) characteristic of de- 
generate stars. 

(ii) On the long-term scale, the outburst characteristics 
change considerably for the massive WD (Model A), but 
only slightly for the low-mass one (Model B). In particular: 

(a) The accreted and ejected masses increase with time, 
as the WD temperature decreases. 

(b) The heavy element abundance (Z) decreases with 
time for Model A, while it oscillates around a constant 
value for model B. 

(c) The helium mass fraction increases with time for 
Model A, while again, it oscillates around a constant value 
for Model B. 

(d) Isotopic ratios settle eventually into constant val- 
ues for Model A and are practically constant throughout 

for Model B. In both cases i^C and "O are overabundant 
as compared to solar values, while is underabundant; 
^^N, however, is overabundant in one case (A) and under- 
abundant in the other (B). 

(iii) The changes in physical characteristics axe not nec- 
essarily monotonic; some, such as the expansion velocity, go 
through an extremum. Consequently, one should be care- 
ful when using parametrised grids of models and interpolat- 
ing between them. Nevertheless, the agreement between the 
long-term calculations and the parametrised models is quite 
good for many of the nova characteristics. 

(iv) Both models settle, eventually, into almost steady 
state where both the WD properties and the outburst char- 
acteristics remain almost unchanged with many repeated 
cycles. However, this state is reached after hundreds of cy- 
cles, which means a long period of time, or a considerable 
change in the mass of the donor star. Although our calcula- 
tions assume a constant accretion rate throughout, it is to 
be expected that spurious dynamical effects on the binary 
orbit that are bound to occur on a long timescale, or the 
change in donor mass, or both, will affect the mass transfer 
rate significantly. A different accretion rate will tend to a 
different steady state. It is thus possible that, in reality, the 
system will never be able to achieve steady state. 

(v) For the first time in numerical modelling of nova out- 
bursts, it is possible to estimate the accuracy (or reliability) 
of the results. The numerical "noise" that appears in the 
various plots may be statistically interpreted as error bars. 
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Of the 3 independent nova parameters, only 2 are left 
— the WD mass and the accretion rate — while the third — 
the WD temperature — is determined by the evolutionary 
course of the accreting WD in the nova system. The present 
calculations were based on the arbitrary assumption that the 
accretion rate is constant (and prescribed). This assumption 
is not realistic; the evolutionary course of the nova system, 
taking into account the interaction between its components, 
determines the mass transfer rate and its evolution, and fu- 
ture studies should follow the evolution of both members 
of the binary system consistently. This leaves the WD mass 
as the single truly independent parameter of nova outburst 
evolution. 
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